// eee.Sheffield.PZ.Math
//
// Copyright ?Ping Zou, 2007
// sg71.cherub@gmail.com

using System;
using System.Collections.Generic;
using System.Text;
using System.IO;
using System.Runtime.Serialization.Formatters.Binary;
using System.Runtime.Serialization;

using eee.Sheffield.PZ.Math;

namespace eee.Sheffield.PZ.Imaging
{
    /// <summary>
    /// GVF Thining sampler, RJMCMC
    /// </summary>
    public class MCMCGVFThinningSampler
    {
        #region Fields
        private StreamWriter _logWriter = null;
        
        private bool _isContinue = false;
        private string _mcmcStateFileName = "";
        
        private MCMCState _mcmcState = null;

        // basic RJMCMC 
        //private TargetDistribution _GMTargetDistribution;
        private LineSegmentConfiguration _lineSegmentListPool;  // line segment pool, data term should have been assigned.               
        //private LineSegmentConfiguration _sourceLineSegmentList;    // source line segment, equals to full configuration        
        private int _totalLineSegment;               // total # of line segments

        // RJMCMC parameters
        private int _total;     // number of MCMC iterations        
        private int _burnin;    // number of burn-in      
        private int _space;     // sample rate
        //private List<LineSegmentConfiguration> _targetSamples;  // target samples, after burn-in, sampled at space
        //private List<LineSegmentConfiguration> _totalSamples; // all samples, including burn-in

        // RJMCMC performance variable
        //private double _sumAcceptanceRate;
        //private double _averageAcceptanceRate;
        //private int _hits;

        // results
        //private PZMath_matrix _lineSegmentHittingProbabilityMap;            // line segments' hitting probability map, gray-level
        //private LineSegmentConfiguration _lineSegmentConfigurationWithMinU; // line segment configuration with max H()        
        //private PZMath_vector _lineSegmentNumberArray;   // number of line segments v.s. MCMC iteration

        // debug variables
        //private double _sumn = 0.0;
        //private double _averagen = 0.0;
        //private double _sumnf = 0.0;
        //private double _averagenf = 0.0;
        //private double _sumnfe = 0.0;
        //private double _averagenfe = 0.0;
        //private double _sumpC = 0.0;
        //private double _averagepC = 0.0;
        //private double _sumpL = 0.0;
        //private double _averagepL = 0.0;
        //private double _sumpG = 0.0;
        //private double _averagepG = 0.0;
        //private double _sumpI = 0.0;
        //private double _averagepI = 0.0;
        //private double _sumhp = 0.0;
        //private double _averagehp = 0.0;
        //private double _sumhd = 0.0;
        //private double _averagehd = 0.0;
        //private double _sumh = 0.0;
        //private double _averageh = 0.0;
        //private int _indexMax = 0;
        #endregion

        #region Properties
        public StreamWriter LogWriter { get { return _logWriter; } set { _logWriter = value; } }

        public bool IsCoutinue { set { _isContinue = value; } }
        public string MCMCStateFileName {set {_mcmcStateFileName = value;}}
        public MCMCState MCMCStateInstance { get { return _mcmcState; } }

        //public TargetDistribution GMTargetDistribution { get { return _GMTargetDistribution; } set { _GMTargetDistribution = value; } }
        public LineSegmentConfiguration LineSegmentListPool { get { return _mcmcState.Pool; } set { _mcmcState.Pool = value; } }
        //public LineSegmentConfiguration SourceLineSegmentList { get { return _sourceLineSegmentList; } set { _sourceLineSegmentList = value; } }
        public int TotalLineSegment { get { return _totalLineSegment; } set { _totalLineSegment = value; } }
        
        // RJMCMC parameters
        public int Total { get { return _total; } set { _total = value; } }
        public int BurnIn { get { return _burnin; } set { _burnin = value; } }
        public int Space { get { return _space; } set { _space = value; } }
        //public List<LineSegmentConfiguration> TargetSamples { get { return _targetSamples; } set { _targetSamples = value; } }
        //public List<LineSegmentConfiguration> TotalSamples { get { return _totalSamples; } set { _totalSamples = value; } }
        
        // RJMCMC performance variables
        //public double SumAcceptanceRate { get { return _sumAcceptanceRate; } set { _sumAcceptanceRate = value; } }
        //public double AverageAcceptanceRate { get { return _averageAcceptanceRate; } set { _averageAcceptanceRate = value; } }
        //public int Hits { get { return _hits; } set { _hits = value; } }
        
        // results
        public List<int> LineSegmentNumberArray { get { return _mcmcState.LineSegmentNumberList; } }
        //public PZMath_matrix LineSegmentHittingProbabilityMap { get { return _lineSegmentHittingProbabilityMap; } }
        public LineSegmentConfiguration LineSegmentConfigurationWithMinU { get { return _mcmcState.UMinLineSegmentConfiguration; } }
        public LineSegmentConfiguration HitRecord { get { return _mcmcState.HitRecord; } }

        // debug variables        
        //public double Sumn { get { return _sumn; } set { _sumn = value; } }
        //public double Averagen { get { return _averagen; } set { _averagen = value; } }
        //public double Sumnf { get { return _sumnf; } set { _sumnf = value; } }
        //public double Averagenf { get { return _averagenf; } set { _averagenf = value; } }
        //public double Sumnfe { get { return _sumnfe; } set { _sumnfe = value; } }
        //public double Averagenfe { get { return _averagenfe; } set { _averagenfe = value; } }
        //public double SumpC { get { return _sumpC; } set { _sumpC = value; } }
        //public double AveragepC { get { return _averagepC; } set { _averagepC = value; } }
        //public double SumpL { get { return _sumpL; } set { _sumpL = value; } }
        //public double AveragepL { get { return _averagepL; } set { _averagepL = value; } }
        //public double SumpG { get { return _sumpG; } set { _sumpG = value; } }
        //public double AveragepG { get { return _averagepG; } set { _averagepG = value; } }
        //public double SumpI { get { return _sumpI; } set { _sumpI = value; } }
        //public double AveragepI { get { return _averagepI; } set { _averagepI = value; } }
        //public double Sumhp { get { return _sumhp; } set { _sumhp = value; } }
        //public double Averagehp { get { return _averagehp; } set { _averagehp = value; } }
        //public double Sumhd { get { return _sumhd; } set { _sumhd = value; } }
        //public double Averagehd { get { return _averagehd; } set { _averagehd = value; } }
        //public double Sumh { get { return _sumh; } set { _sumh = value; } }
        //public double Averageh { get { return _averageh; } set { _averageh = value; } }
        //public int IndexMax { get { return _indexMax; } }
        #endregion

        #region Constructor
        /// <summary>
        /// number of iterations, number of burn-in, sample rate
        /// </summary>
        //public MCMCGVFThinningSampler(int total, int burnin, double lamda, double w, double h)
        public MCMCGVFThinningSampler(int total, int burnin, int space, LineSegmentConfiguration pool)
        {
            _total = total;
            _burnin = burnin;
            _space = space;

            if (_burnin >= _total)
                throw new ApplicationException("GVFThinningSampler::GVFThinningSampler(), burn-in length is great than sample length!");            
            
            // line segment pool            
            _lineSegmentListPool = new LineSegmentConfiguration(pool); // new LineSegmentConfiguration(pool);                        
            _totalLineSegment = pool.N;
        } // GVFThinningSampler()
        #endregion

        #region RJMCMC sampler GeyerMoller algorithm        

        /// <summary>
        /// U(), H() = exp(-U())
        /// </summary>
        /// <returns></returns>
        public static double U(LineSegmentConfiguration x)
        {
            bool slow = true;
            // calculate hp
            x.FreeEndsFreeSegments();
            if (slow)
                x.JoinSegments();     // undo join to save computation time
            x.CalculatePC();
            x.CalculateSumPC();

            if (slow)
            {
                x.CalculateTotalL();
                x.CalculateAverageL();
            }

            // explicitly calculate sum(pL(s)), sum(pG(s)), and sum(pI(s)) for the efficiency sake
            int count = x.Configure.Count;
            double sumPL = 0.0;
            double sumPG = 0.0;
            double sumPI = 0.0;
            for (int s = 0; s < count; s++)
            {
                if (slow)
                {
                    // calculate sum(pL(s))
                    x.Configure[s].CalculatePL(x.TotalLength, x.AverageLength);

                    // try different weight
                    double pL = x.Configure[s].PL;
                    if (pL < 0.0)
                        sumPL += pL * 100;
                    else
                        sumPL += pL;
                    x.SumPL = sumPL;
                }

                // calculate sum(pG(s))
                x.Configure[s].CalculatePG(x.Gt);
                double pG = x.Configure[s].PG;
                if (pG < 0.0)
                    sumPG += pG * 100;
                else
                    sumPG += pG;
                x.SumPG = sumPG;

                // calculate sum(pI(s))
                x.Configure[s].CalculatePI(x.It, x.KI);
                double pI = x.Configure[s].PI;
                if (pI < 0.0)
                    sumPI += pI * 10000;
                else
                    sumPI += pI;
                x.SumPI = sumPI;
            }

            double up = 0;
            if (slow)
                up = x.W1 * x.N + x.W2 * x.NF + x.W3 * x.NFE + x.WC * x.SumPC +x.WL * x.SumPL;
            else
                up = x.W1 * x.N + x.W2 * x.NF + x.W3 * x.NFE + x.WC * x.SumPC;
            
            // calculate hd           
            double ud = x.WG * x.SumPG + x.WI * x.SumPI;

            // h = hp + hd;
            double u = up + ud;           
                        
            return u;
        } // U()        

        /// <summary>
        /// birth accept probability (birth accept rate)
        /// R = h(y) * (N - n(x)) / h(x) / n(y)
        /// log(R) = -u(y) + log((N - n(x)) - (-u(x)) - log(n(y))
        /// </summary>
        /// <returns></returns>
        private double BirthAcceptProbability(LineSegmentConfiguration x, LineSegmentConfiguration y, out double UXn)
        {
            // calculate birth accept probability 
            double R;
            double logR;
            double uY = -1.0 * U(y);
            double uX = -1.0 * U(x);
            UXn = uX;
            //R = hY * (_totalLineSegment - x.N) / hX / y.N;
            logR = uY + System.Math.Log(_totalLineSegment - x.N) - uX - System.Math.Log(y.N);
            
            if (logR > 0)
                R = 1.0;
            else
                R = System.Math.Exp(logR);
            return R;
        } // BirthAcceptProbability()

        /// <summary>
        /// death accept probability
        /// R = h(y) * n(x) / h(x) / (N - n(y))
        /// logR = u(y) + log(n(x)) - u(x) - log(N - n(y))
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <returns></returns>
        private double DeathAcceptProbability(LineSegmentConfiguration x, LineSegmentConfiguration y, out double UXn)
        {
            double R;
            double logR;

            double uY = -1.0 * U(y);
            double uX = -1.0 * U(x);
            UXn = uX;

            //R = hY * x.N / hX / (_totalLineSegment - y.N);
            logR = uY + System.Math.Log(x.N) - uX - System.Math.Log(_totalLineSegment - y.N);
            
            if (logR > 0)
                R = 1.0;
            else
                R = System.Math.Exp(logR);
            return R;
        } // DeathAcceptProbability()

        /// <summary>
        /// do accept rejection sampling
        /// </summary>
        /// <returns></returns>
        public void Sample()
        {
            // RJMCMC starts from full line segment configuration            
            LineSegmentConfiguration start = new LineSegmentConfiguration(_lineSegmentListPool);
            _lineSegmentListPool.Empty();       // pool now is empty

            // initial work space
            long seed = DateTime.Now.Millisecond;
            ParkMillerUniform uniformRandom = new ParkMillerUniform(seed);

            if (_isContinue && File.Exists(_mcmcStateFileName))
            {
                _mcmcState = (MCMCState)ReadFile(_mcmcStateFileName);
                _total += _mcmcState.N;
            }
            else
            {
                _mcmcState = new MCMCState();
                _mcmcState.HitRecord = new LineSegmentConfiguration(start);
                _mcmcState.Xn = start;
                _mcmcState.Pool = new LineSegmentConfiguration(_lineSegmentListPool);
                _mcmcState.LineSegmentNumberList = new List<int>(_total);
            }

            while (_mcmcState.N < _total)
            {
                int xLength = _mcmcState.Xn.N;
                _mcmcState.IsDeath = false;
                _mcmcState.IsRejected = false;                

                /// Birth or Death ?
                #region Birth or Death ?
                _mcmcState.BirthProb = uniformRandom.Sample();
                if (xLength > 0 && xLength < _totalLineSegment)
                    _mcmcState.DeathProb = 1 - _mcmcState.BirthProb;
                else if (xLength == 0)
                // only birth
                {
                    _mcmcState.BirthProb = 1.0;
                    _mcmcState.DeathProb = 0.0;
                }
                else
                // only death
                {
                    _mcmcState.DeathProb = 1.0;
                    _mcmcState.BirthProb = 0.0;
                }
                #endregion

                double uXn;
                if (_mcmcState.BirthProb > _mcmcState.DeathProb)
                // Birth
                {
                    _mcmcState.IsDeath = false;
                    // sample a new line segment from pool
                    _mcmcState.PoolIndex = (int)System.Math.Floor(uniformRandom.Sample() * _mcmcState.Pool.N);
                    LineSegment addedLineSegment = _mcmcState.Pool.Configure[_mcmcState.PoolIndex];

                    // propose Y                   
                    _mcmcState.Y = new LineSegmentConfiguration(_mcmcState.Xn);
                    _mcmcState.Y.Add(addedLineSegment);

                    _mcmcState.AcceptRate = BirthAcceptProbability(_mcmcState.Xn, _mcmcState.Y, out uXn);
                    _mcmcState.SumAcceptRate += _mcmcState.AcceptRate;
                    double u = uniformRandom.Sample();
                    if (u <= _mcmcState.AcceptRate)
                    // accept Y
                    {
                        _mcmcState.IsRejected = false;
                        _mcmcState.StayCount = 0;

                        _mcmcState.XnA1 = new LineSegmentConfiguration(_mcmcState.Y);

                        if (_mcmcState.N >= _burnin)
                        {
                            if (_mcmcState.SpaceCount >= _space)
                            {
                                // count hits of line segments
                                int sourceLength = _mcmcState.HitRecord.N;
                                for (int s = 0; s < xLength; s++)
                                {
                                    int lineSegmentIndex = _mcmcState.Xn.Configure[s].Index;
                                    _mcmcState.HitRecord.Configure[lineSegmentIndex].Hit++;
                                }
                                _mcmcState.SpaceCount = 0;
                            }
                            else
                                _mcmcState.SpaceCount++;
                        }
                        // remove l from pool
                        _mcmcState.Pool.RemoveAt(_mcmcState.PoolIndex);
                    }
                    else
                    // reject Y
                    {
                        _mcmcState.IsRejected = true;
                        _mcmcState.StayCount++;

                        _mcmcState.XnA1 = new LineSegmentConfiguration(_mcmcState.Xn);

                        if (_mcmcState.N >= _burnin)
                        {
                            if (_mcmcState.SpaceCount >= _space)
                            {
                                // count hits of line segments
                                int sourceLength = _mcmcState.HitRecord.N;
                                for (int s = 0; s < xLength; s++)
                                {
                                    int lineSegmentIndex = _mcmcState.Xn.Configure[s].Index;
                                    _mcmcState.HitRecord.Configure[lineSegmentIndex].Hit++;
                                }
                                _mcmcState.SpaceCount = 0;
                            }
                            else
                                _mcmcState.SpaceCount++;
                        }
                    }
                }
                else
                // Death
                {
                    _mcmcState.IsDeath = true;
                    // uniformly remove v from X
                    _mcmcState.RemovedIndex = (int)System.Math.Floor(uniformRandom.Sample() * xLength);
                    // propose Y                    
                    _mcmcState.Y = new LineSegmentConfiguration(_mcmcState.Xn);

                    LineSegment removedLineSegment = _mcmcState.Y.Configure[_mcmcState.RemovedIndex];
                    _mcmcState.Y.RemoveAt(_mcmcState.RemovedIndex);

                    // accept probability
                    _mcmcState.AcceptRate = DeathAcceptProbability(_mcmcState.Xn, _mcmcState.Y, out uXn);
                    _mcmcState.SumAcceptRate += _mcmcState.AcceptRate;
                    double u = uniformRandom.Sample();
                    if (u <= _mcmcState.AcceptRate)
                    // accept Y
                    {
                        _mcmcState.IsRejected = false;
                        _mcmcState.StayCount = 0;

                        _mcmcState.XnA1 = new LineSegmentConfiguration(_mcmcState.Y);

                        if (_mcmcState.N >= _burnin)
                        {
                            if (_mcmcState.SpaceCount >= _space)
                            {
                                _mcmcState.SpaceCount = 0;

                                // count hits of line segments
                                int sourceLength = _mcmcState.HitRecord.N;
                                for (int s = 0; s < xLength; s++)
                                {
                                    int lineSegmentIndex = _mcmcState.Xn.Configure[s].Index;
                                    _mcmcState.HitRecord.Configure[lineSegmentIndex].Hit++;
                                }

                            }
                            else
                                _mcmcState.SpaceCount++;
                        }
                        
                        // add l back to pool
                        _mcmcState.Pool.Add(removedLineSegment);
                    }
                    else
                    // reject Y
                    {
                        _mcmcState.IsRejected = true;
                        _mcmcState.StayCount++;

                        _mcmcState.XnA1 = new LineSegmentConfiguration(_mcmcState.Xn);

                        //if (n < Nm1)
                        //    _totalSamples.Add(X);

                        if (_mcmcState.N >= _burnin)
                        {
                            if (_mcmcState.SpaceCount >= _space)
                            {
                                _mcmcState.SpaceCount = 0;

                                // count hits of line segments
                                int sourceLength = _mcmcState.HitRecord.N;
                                for (int s = 0; s < xLength; s++)
                                {
                                    int lineSegmentIndex = _mcmcState.Xn.Configure[s].Index;
                                    _mcmcState.HitRecord.Configure[lineSegmentIndex].Hit++;
                                }
                            }
                            else
                                _mcmcState.SpaceCount++;
                        }                        
                    }
                }                

                // shuffle if stay count over 200
                if (_mcmcState.StayCount > 100)
                {
                    _mcmcState.StayCount = 0;

                    // shufffle the configuration
                    // randomly choose one line segment from the pool, replace the random one in the current configuration
                    // do the above n(x) times
                    LineSegmentConfiguration currentLSC = _mcmcState.XnA1;
                    int currentLineSegments = currentLSC.N;

                    for (int i = 0; i < currentLineSegments; i++)
                    {
                        // randomly choose one line segment from the pool
                        int poolIndex = (int)System.Math.Floor(uniformRandom.Sample() * MCMCStateInstance.Pool.N);
                        LineSegment addedLineSegment = _mcmcState.Pool.Configure[poolIndex];

                        // uniformly remove line segment from currentLSC
                        int removeIndex = (int)System.Math.Floor(uniformRandom.Sample() * currentLineSegments);
                        LineSegment removedLineSegment = currentLSC.Configure[removeIndex];

                        // replace
                        currentLSC.RemoveAt(removeIndex);
                        currentLSC.Add(addedLineSegment);

                        // save removed line segment to the pool
                        _mcmcState.Pool.Add(removedLineSegment);
                    }
                    
                    Log("---- Shuffled!!!! ----");
                }
               
                // summary statistics                
                _mcmcState.AverageAcceptRate = _mcmcState.SumAcceptRate / ((double)_mcmcState.N + 1.0);

                double uX = -1.0 * uXn; // U(_mcmcState.Xn);

                StringBuilder tempLogStringBuilder = new StringBuilder();

                if (uX < MCMCStateInstance.UMin)
                {
                    _mcmcState.UMin = uX;
                    _mcmcState.UMinLineSegmentConfiguration = new LineSegmentConfiguration(_mcmcState.Xn);
                    //_indexMax = _mcmcState.N - 1;
                }
                tempLogStringBuilder.Append("n = " + (MCMCStateInstance.N));

                if (MCMCStateInstance.IsDeath)
                    tempLogStringBuilder.Append(" Death");
                else
                    tempLogStringBuilder.Append(" Birth");
                if (MCMCStateInstance.IsRejected)
                    tempLogStringBuilder.Append(" Reject");
                else
                    tempLogStringBuilder.Append(" Accept");
                tempLogStringBuilder.Append(" # = " + _mcmcState.Xn.N);
                _mcmcState.LineSegmentNumberList.Add(_mcmcState.Xn.N);
                tempLogStringBuilder.Append(" u(x) = " + String.Format("{0:0.0000e+00}", uX));
                tempLogStringBuilder.Append(" accept rate = " + String.Format("{0:0.0000e+00}", _mcmcState.AcceptRate));
                tempLogStringBuilder.Append(" average accept rate " + String.Format("{0:f}", _mcmcState.AverageAcceptRate));
                Log(tempLogStringBuilder.ToString());

                _mcmcState.N++;
                _mcmcState.Xn = _mcmcState.XnA1;

                if (_mcmcStateFileName != null)
                    _mcmcState.WriteFile(_mcmcStateFileName);

            } // end of while

            Log("ux min = " + String.Format("{0:0.0000}", _mcmcState.UMin)); // + " index = " + _indexMax);
            _mcmcState.AverageAcceptRate = _mcmcState.SumAcceptRate / (double)_total;

            //_lineSegmentNumberArray = GetLineSegmentNumberArray(_targetSamples);
            //_lineSegmentConfigurationWithMinU = GetLineSegmentConfigurationWithMinU(_targetSamples);

            //return _targetSamples;
        } // Sample()
        #endregion

        #region Log method
        /// <summary>
        /// Log method
        /// </summary>
        /// <param name="logText"></param>
        private void Log(string logText)
        {
            Console.WriteLine(logText);
            if (_logWriter != null)
            {
                _logWriter.WriteLine(logText);
                _logWriter.Flush();
            }
        } // Log()
        #endregion

        #region IO methods

        #endregion

        #region Summary statistics methods
        /// <summary>
        /// get line segement configuration with max H
        /// </summary>
        /// <param name="lineSegmentationConfigurationSamples"></param>
        /// <returns></returns>
        private LineSegmentConfiguration GetLineSegmentConfigurationWithMinU(List<LineSegmentConfiguration> lineSegmentationConfigurationSamples)
        {
            int sampleCount = lineSegmentationConfigurationSamples.Count;
            double minU = Double.MaxValue;
            int minIndex = 0;
            double u;
            for (int i = 0; i < sampleCount; i++)
            {
                u = U((LineSegmentConfiguration)lineSegmentationConfigurationSamples[i]);
                if (u < minU)
                {
                    minU = u;
                    minIndex = i;
                }
            }

            LineSegmentConfiguration dstLineSegmentConfiguration = new LineSegmentConfiguration((LineSegmentConfiguration)lineSegmentationConfigurationSamples[minIndex]);
            //dstLineSegmentConfiguration.ClearLoops();
            //dstLineSegmentConfiguration.ClearFreeLineSegments();
            return dstLineSegmentConfiguration;
        } // GetLineSegmentConfigurationWithMinU()

        //public void GetLineSegmentConfigurationWithMinU()
        //{
        //    _lineSegmentConfigurationWithMinU = GetLineSegmentConfigurationWithMinU(_targetSamples);            
        //} // GetLineSegmentConfigurationWithMaxH()

        /// <summary>
        /// get line segment number array v.s. MCMC iteration
        /// </summary>
        /// <param name="lineSegmentConfigurationSamples"></param>
        /// <returns></returns>
        private PZMath_vector GetLineSegmentNumberArray(List<LineSegmentConfiguration> lineSegmentConfigurationSamples)
        {
            int sampleCount = lineSegmentConfigurationSamples.Count;
            PZMath_vector lineSegmentNumberArray = new PZMath_vector(sampleCount);
            for (int i = 0; i < sampleCount; i++)
            {
                lineSegmentNumberArray[i] = ((LineSegmentConfiguration)lineSegmentConfigurationSamples[i]).N;
            }

            return lineSegmentNumberArray;
        } // GetLineSegmentNumberArray()

        //public void GetLineSegmentNumberArray()
        //{
        //    _lineSegmentNumberArray = GetLineSegmentNumberArray(_targetSamples);
        //} // GetLineSegmentNumberArray()
        #endregion

        #region IO methods
        /// <summary>
        /// load MCMCState from file
        /// </summary>
        /// <param name="fileName"></param>
        /// <returns></returns>
        private MCMCState ReadFile(string fileName)
        {
            Stream streamRead = File.OpenRead(fileName);
            BinaryFormatter binaryRead = new BinaryFormatter();
            MCMCState mcmcState = (MCMCState)binaryRead.Deserialize(streamRead);
            streamRead.Close();
            return mcmcState;
        } // ReadFile()
        #endregion

        #region Example
        #endregion
    }

    /// <summary>
    /// MCMC state class
    /// </summary>
    [Serializable]
    public class MCMCState
    {
        #region Fields
        private LineSegmentConfiguration _Xn = null;
        private LineSegmentConfiguration _XnA1 = null;
        private LineSegmentConfiguration _Y = null;
        private LineSegmentConfiguration _hitRecord = null;
        private LineSegmentConfiguration _pool = null;
        private LineSegmentConfiguration _uMinLineSegmentConfiguration = null;

        private int _n;                 // MCMC iteration
        private int _spaceCount;        // sample space count
        private int _stayCount;         // stay count, for shuffle purpose
        private double _uMin = Double.MaxValue;           // min U(X)

        private double _birthProb;      // birth prob
        private double _deathProb;      // death prob = 1 - mcmcState.BirthProb
        private int _poolIndex;      // random index in pool
        private int _removedIndex;    // random removed index
        private double _acceptanceRate; // accept rate
        private double _averageAcceptanceRate;  //average accept rate
        private double _sumAcceptanceRate;      // sum accept rate

        private bool _isDeath;          // is death/birth?
        private bool _isRejected;       // is Y accepted/rejected?

        private List<int> _lineSegmentNumberList;   // number of line segment v.s. MCMC iteration
        #endregion

        #region Properties
        public LineSegmentConfiguration Xn { get { return _Xn; } set { _Xn = value; } }
        public LineSegmentConfiguration XnA1 { get { return _XnA1; } set { _XnA1 = value; } }
        public LineSegmentConfiguration Y { get { return _Y; } set { _Y = value; } }
        public LineSegmentConfiguration Pool { get { return _pool; } set { _pool = value; } }
        public LineSegmentConfiguration UMinLineSegmentConfiguration { get { return _uMinLineSegmentConfiguration; } set { _uMinLineSegmentConfiguration = value; } }
        public LineSegmentConfiguration HitRecord { get { return _hitRecord; } set { _hitRecord = value; } }

        public int N { get { return _n; } set { _n = value; } }
        public int SpaceCount { get { return _spaceCount; } set { _spaceCount = value; } }
        public int StayCount { get { return _stayCount; } set { _stayCount = value; } }
        public double UMin { get { return _uMin; } set { _uMin = value; } }

        public double BirthProb { get { return _birthProb; } set { _birthProb = value; } }
        public double DeathProb { get { return _deathProb; } set { _deathProb = value; } }
        public int PoolIndex { get { return _poolIndex; } set { _poolIndex = value; } }
        public int RemovedIndex { get { return _removedIndex; } set { _removedIndex = value; } }
        public double AcceptRate { get { return _acceptanceRate; } set { _acceptanceRate = value; } }
        public double AverageAcceptRate { get { return _averageAcceptanceRate; } set { _averageAcceptanceRate = value; } }
        public double SumAcceptRate { get { return _sumAcceptanceRate; } set { _sumAcceptanceRate = value; } }

        public bool IsDeath { get { return _isDeath; } set { _isDeath = value; } }
        public bool IsRejected { get { return _isRejected; } set { _isRejected = value; } }

        public List<int> LineSegmentNumberList { get { return _lineSegmentNumberList; } set { _lineSegmentNumberList = value; } }
        #endregion

        /// <summary>
        /// save MCMCState to file
        /// </summary>
        /// <param name="fileName"></param>
        /// <param name="mcmcState"></param>
        public void WriteFile(string fileName)
        {
            Stream streamWrite = File.Create(fileName);
            BinaryFormatter binaryWrite = new BinaryFormatter();
            binaryWrite.Serialize(streamWrite, this);
            streamWrite.Close();
        } // WriteFile
    }
}
